Effects of disorder on atomic density waves and spin-singlet dimers in one-dimensional 

optical lattices 
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Using the Bethe-ansatz density-functional theory, we study a one- dimensional Hubbard model 
of confined attractively interacting fermions in the presence of a uniformly distributed disorder. 
The strongly-correlated Luther-Emery nature of the attractive one-dimensional Hubbard model is 
fully taken into account as the reference system in the density-functional theory. The effects of 
the disorder are investigated on the atomic-density waves in the weak-to-intermediate attractive 
interaction and on the spin-singlet dimers of doubly occupied sites in the strongly attractive regime. 
It is found that atomic density waves are sensitive to the disorder and the spin-singlet dimers of 
doubly occupied sites are quite unstable against the disorder. We also show that very weak disorder 
could smear the singularities in the stiffness, thus suppress the spin-singlet pairs. 
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I. INTRODUCTION 

Disorder plays an important role in physical proper- 
ties of quantum many-body systems, particularly in solid 
state physics. Examples of phenomena related to disor- 
der include a drastic suppression of superfluidity of 4 He 
in nanoporous media,— an apparent metal-insulator quan- 
tum phase transition in a two-dimensional electron liquid 
induced by the interplay between long-range Coulomb in- 
teractions and disorder,- a strong dependence of super- 
solidity on the amount of disorder in solid 4 He^, and a 
halted expansion of a Bose-Einstein condensate (BEC) 
in a random potential. 4 

Intriguing achievements of ultracold atomic systems in 
optical lattices have opened an exciting field in manipu- 
lating many of the model Hamiltonian parameters in a 
clean and controlled way. For example, experimentalists 
can control the on-site interaction strength, confining po- 
tentials, the effective dimensionality, and different forms 
of disorder to simulate phenomena in solid-state physics 
like the superfluid-Mott-insulator transitions and BEC- 
BCS crossover of long range phase coherence of fermionic 
pairs. 6 The on-site interaction can be tuned either indi- 
rectly by changing the strength of the lasers that create 
the optical lattice potential or directly by using a Fes- 
hbach resonance, which allows to change the scattering 
length from — oo to oo. Optical lattices can be used to 
find new exotic quantum states of matter, for example, 
by trapping Bose-Fermi mixtures^ and multiple hyper- 
fine states^ that have no direct analogues in condensed 
matter systems. Besides the clean realizations of many 
condensed-matter lattice models in optical lattices, dis- 
order can also be generated in optical lattices. 

Disorder in optical lattices is attracting increasing at- 
tention of many theorists and experimentalists due to 
its unique flexibility and unprecedented controllability 
available in producing and observing disordered quantum 
degenerate gases. The unique controllability gains the 
research on optical lattices of disorder much advantage 
over the disorder research in electronic systems, where 



the disorder is fixed by a specific realization of the sam- 
ple and the interaction is deemed to be the long-range 
Coulomb one3 Disorder in optical lattices can be cre- 
ated in different ways: an unbounded speckle pattern 
produced by shining a laser beam through a diffusive 
plate,— ii£ bounded incommensurate bichromatic lattices 
produced by combining the primary optical lattice with 
a secondary lattice (not random but completely deter- 
ministic quasidisorder ) 1 1 1 1 1 2 1 1 3 and a superimposed dis- 
order*^ Other novel kinds of disorder have been pro- 
posed, for example, a bounded disorder in the strength 
of interatomic interactions, which can be realized near 
a Feshbach resonance^ and randomly located impurity 
atoms with a different internal state trapped at the nodes 
of an optical latticed Several proposals have been put 
forward in regard to disordered ultracold bosonic atom 
gases, where a rich variety of exotic new quantum phases 
like the Bose glass, the Bose-Anderson glass, and the 
Mott glass have been predicted i 14 ' 17 i 18 Anderson local- 
ization is predicted to be observable from the expansion 
of an initially confined interacting one-dimensional (ID) 
BEC in a weak disorder^ 

Cold Bose gases have been successfully trapped in ID 
geometries. A Tonks- Girardeau gas of bosonic 87 Rb 
atoms 2 ^ was realized experimentally in a ID optical lat- 
tice. A 87 Rb atomic gas inside a disordered ID optical 
lattice (OL) has been used to study the interplay be- 
tween repulsive interactions and disorder^ Towards the 
experimental test of predicted exotic phases, the Florence 
group gave us the first experimental hint of Bose glass 
by measuring the excitation energy spectrum and coher- 
ence i^ 3 - In the case of fermionic atom gases, however, it is 
more difficult to cool them down because of their Pauli- 
limited collision rate, different from the s-wave collisional 
pattern of bosonic atoms. Up to now a two-component 
Fermi gas of K atoms has been prepared in a quasi- 
1D geometryjSi An interacting Fermi gas of 40 K atoms 
has also been demonstrated in three-dimensional optical 
lattices^ We think that trapping interacting gases of 
fermionic atoms in ID optical lattices is within the reach 
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of the present-day techniques. 

These recent developments in experiment make an ex- 
ploration of fermions in disordered optical lattices timely. 
Some proposals related to fermionic atom gases in dis- 
ordered OLs have already been made to revisit disor- 
dered condensed matter models. For example, in the 
work of Paredes et al.^ it was reported that fermionic 
atoms in optical superlattices exhibit strongly correlated 
phenomena, from Kondo singlet formation to magnetism 
of localized spins. Yamashita et al. studied the char- 
acteristic properties of fermions trapped in a ID opti- 
cal superlattice with 2-site periodicity and found that 
three different insulating phases (of band-, bond-charge- 
density-wave- and Mott-type) can emerge.— Repulsive 
interacting Fermi gases in disordered ID OLs have been 
studied in Ref. 2j| where the effects of disorder on local 
Mott-insulating and band-insulating regions have been 
analyzed in detail. 

Motivated by these experimental and theoretical sce- 
nario, in this paper we study the interplay between at- 
tractive interactions and uniform randomness in a two- 
component Fermi gas loaded in a ID OL under har- 
monic confinement. We investigate systematically the 
effects of a random potential on these systems making 
use of Bethe-ansatz density-functional theory . 26 i 27 In the 
clean limit, the system exhibits atomic-density waves for 
weak-to-intermediate attractive interactions and band- 
insulating regions for strongly attractive interactions^ 

The structure of the paper is the following. In Sec. II, 
the model Hamiltonian is briefly introduced. In Sec. Ill, 
we present a lattice version of density-functional theory 
based on the exactly solvable ID homogeneous Hubbard 
model. Finally, numerical results and some conclusions 
are shown in Sec. IV. 



II. THE ID RANDOM ATTRACTIVE 
FERMI-HUBBARD MODEL 

A ID Fermi gas trapped by an external harmonic po- 
tential and by a stationary optical potential can be de- 
scribed by a ID single-band Fcrmi-Hubbard model, if wc 
assume that the lattice potential is deep enough and that 
the energy separation between the first and the second 
band of the lattice is much larger than other energy scales 
involved: 

H s = -t (c\ a c l+la + H.c.) + hi-\hii 

i—l,(7 i—1 

+V 2 Y j {i-N s /2fn i . (1) 

i=l 

Here a =T,I is a pseudospin-1/2 label for two internal 

hyperfine states, hi — ™i<? = So- ^la^ is t ne total 
site occupation operator, t is the tunneling between near- 
est neighbors, U is the on-site attractive interaction and 
V2 is the strength of harmonic potential. A system of 



N a lattice sites and Nt interacting fermions is consid- 
ered. The ID random attractive Fermi-Hubbard model 
is written as, 



JV, 

H = H s +Y^e l h l . (2) 

i=i 



The effect of disorder is simulated by the last term in 
Eq. @, where is randomly chosen at each site with a 
uniform distribution in the range [— W/2,W/2]. We as- 
sume that the lattice is deep enough so that the disorder 
alters only Si but not the strength of tunneling between 
neighboring sites t and the on-site interatomic attractive 
interactions U^L We emphasize here that the attractive 
contact interaction has recently been realized in exper- 
iments with quantum gases by means of Feshbach res- 
onances. 29 Experimentally the uniform random disorder 
can be simulated approximately by using a quasiperiodic 
disorder.— 

Theoretical studies of disordered systems suffer a few 
difficulties. Firstly one has to simulate either small sam- 
ples with a large amount of realizations of disorder, or to 
simulate very large samples. The required self-averaging 
makes sure that the system characterization is indepen- 
dent of one particular disorder realization. Due to the fi- 
nite size of the real OLs, we concentrate in this paper on 
one-dimensional chains of few hundred lattice sites. Ex- 
tension to 1000 lattice sites does not increase the numer- 
ical effort substantially. Secondly the interplay between 
strong interactions and disorder is a true challenge. To 
overcome these difficulties, we apply the Bethe-ansatz- 
based density-functional theory within the local density 
approximation (BALDA), which uses the exactly solv- 
able one-dimensional Luther-Emery liquid as the refer- 
ence system. Computationally, BALDA takes a few sec- 
onds to several minutes instead of a few hours to a few 
days when using the quantum Monte Carlo (QMC) and 
density-matrix renormalization group (DMRG) for a sin- 
gle density profile calculation. But unlike BALDA, QMC 
and DMRG also provide an access to correlation func- 
tions and to the momentum distribution. In this study, 
the ground state properties in the presence of disorder 
are obtained by means of a disorder ensemble average, 
((C))dis = lmiM^oo £) a=1 0(a). For example, the 
site occupation A/i is calculated as, A/j = {{ni))dis, where 
m is the ground-state site occupation defined in Eq. (|3j). 
The highly efficient computational methods allow us to 
perform a large amount of disorder realizations. Aver- 
aging over a large number of realizations is necessary as 
the strength of disorder increases. In this study, we take 
A4 = 10 4 . We find that the density profiles are stable 
against a further increase of disorder realizations. 
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III. LATTICE DENSITY-FUNCTIONAL 
THEORY AND THOMAS-FERMI 
APPROXIMATION 

1. Lattice density-functional theory 

The so-called site-occupation functional theory 
(SOFT)2&2L22 is a powerful tool to calculate the 
ground-state (GS) properties of an inhomogeneous 
lattice Hamiltonian. In our case the inhomogeneity is 
caused by harmonic potential and disorder. Within 
SOFT the exact GS site occupation, 



= 2 



occ. 

El (<*) 
m 



(3) 



can be obtained by solving self-consistently the lattice 
Kohn-Sham (KS) equations 



attractive Hubbard model under harmonic confinement, 
the agreement between the BALDA and DMRG is excel- 
lent in certain range of parameters.— For large u and/or 
small values of V^/t, the BALDA scheme overestimates 
the amplitude of the atomic density waves (ADWs) and 
thus improved exchange-correlation functionals are de- 
manded^ In this paper, we choose interaction strength 
and confining potential within the range of values where 
BALDA works well compared to DMRG. In the strong 
attractive regime of Eq. (Q]), a flat region of doubly- 
occupied sites emerges at the trap center, manifesting 
a state of tightly bound spin-singlet dimers. In this case, 
we resort to Thomas-Fermi approximation, which gives 
an overall shape of the density site occupation and almost 
the same flat region of doubly-occupied sites emerging at 
the trap center but misses the density oscillation at the 
edges and atom tunneling beyond the Thomas-Fermi ra- 
dius. 



E 

3=1 



1,3 



(4) 



Here, the effective KS potential is given by vf s = Un,i/2+ 
vf c + V2(i — N s /2) 2 + Si. The sum over a runs over 
all the occupied orbitals, and factor 2 in Eq. is due 
to spin degeneracies. The first term in the Kohn-Sham 
potential is the Hartree mean-field contribution, while 
t>* c = SE xc [n]/Sni\ GS is the exchange-correlation (xc) po- 
tential, defined as the derivative of the xc energy E KC \n\ 
evaluated at the GS site occupation [the interested reader 
is encouraged to consult Ref. 27, Appendix A for a deriva- 
tion of Eq. flU)]. 

The total GS energy of the system is given by 



E[n] 



«rm-X;tfn?/4- 



-E xc [n}. (5) 



2. Thomas- Fermi approximation 

In the strong attractive limit, we calculate the GS 
properties of H in Eq. (J5]) within LDA through calculat- 
ing the chemical potential of the system from the local 
equilibrium condition, 



M = (J>hom{n,u)\ n ^ n . + V 2 (i - N s /2f + e, 



(8) 



derived from the direct minimization of total energy func- 
tionals, 



E[n] = ^e GS (n,w) 

% n — >rii 



(9) 



In the actual calculation, £ xc [n] has to be approximated. 
In this work we employ a Bethe-ansatz-based local den- 
sity approximation for the xc potential, 



horn 



I BALDA 



(n, u) 



(6) 



where, in analogy with ab initio DFT, the xc potential 
(n, u) of the ID homogeneous Hubbard model is de- 



,hom 



fined by 

i£r(n,«) 



d_ 

On 



,(n,0) 



U 



(7) 



Here, e G s( n j u ) is the GS energy per site of the ID homo- 
geneous system as a function of filling n = Nf/N s and 
interaction strength u = —\U\/t. Thus, within the local 
density approximation scheme proposed in Eq. (|6|), the 
only necessary input is the xc potential of ID homoge- 
neous Hubbard model, which can be found numerically 
from the Bcthe-ansatz equations.— BALDA has been 
shown to provide a very good account of the GS prop- 
erties of ID inhomogeneous lattice systems! 27 ' 28 For the 



The particle number is determined by the normalization 
condition Nf = n i- Here, /Xhom(^) u) is the chemi- 

cal potential of the homogeneous system calculated from 
an appropriate set of Bethe-ansatz equations. This ap- 
proach is termed as the Thomas-Fermi approximation 
(TFA), 27 equivalent to using an LDA also for the nonin- 
teracting kinetic energy functional in the DFT scheme. 
The TFA takes a non-interacting approximation for the 
kinetic energy but considers the exchange-correlation en- 
ergy of the corresponding ID homogeneous system, which 
is incorporated in Hh om (n, u). The performance of TFA 
is compared with the DMRG in the next section for the 
case of large |u|. 



IV. NUMERICAL RESULTS AND DISCUSSION 

Without harmonic confinement and disorder (V2 = 0, 
W = 0), the model described by H has been exactly 
solved by Lieb and Wu using the Bethe-ansatz. For 
attractive interactions it belongs to the Luther-Emery 
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FIG. 1: (color online) Site occupation M% as a function of site 
position i for a system of Nf = 30, N s = 100 and u = —2, 
Vijt — 1.0 x 10 -3 . The inset shows the amplitude of the 
ADW, Aadw, as a function of the strength of the disorder 
W/t for different interaction strengths of u — —2 and u = — 1. 
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FIG. 2: (color online) Fraunhofer structure factor S(q) for 
different disorder strengths in a system of Nf — 30, N 3 = 
100, u — —2, and in the presence of a harmonic potential 
V 2 /t = 1.0 x 10" 3 . 



liquid universality class. Luthcr-Emcry liquids arc Lut- 
tinger liquids exhibiting gapless density excitations and 
gapped spin excitations. Without disorder (W = 0), this 
system exhibits for weak to intermediate attractive inter- 
action compound phases characterized by the coexistence 
of spin pairing and atomic-density waves. — The ap- 
pearance of atomic-density waves in ID Fermi gases with 
attractive interactions under confinement is due to the 
in-phasc Friedcl oscillations of spin-resolved ground-state 
density profiles. For strong atom-atom attractive inter- 
actions, a state of tightly bound spin-singlet dimers of 
doubly occupied sites appears at the center of the trapi 2 ^ 

Here we focus on the interplay between attractive in- 
teractions and disorder. For the weak-to-intermediate at- 
tractive interaction, we perform full DFT calculations for 
disorder realizations up to 10 4 disorder ensemble average. 
In Fig. [TJ we show the GS site occupation for Nf = 30 
atoms in a lattice of N s = 100 sites, inside a trap with 
Vz/t = 10~ 3 . Without disorder, the atomic density wave 
is clearly seen as a consequence of Luther- Emery pairing 
and the formation of stable spin-singlet dimers between 
different pseudospins. For W = 0, the density profile cal- 
culated from BALDA has been compared with the one 
from DMRG and both are in a very good agreement^ At 
W — 0.5 the ADW is largely suppressed by the disorder. 
For large values of \u\, the amplitude of ADWs does not 
scale to zero when spin-singlet dimers occupy the deep- 
est valleys in the disorder landscape due to the attractive 
nature of interaction. As a result, there is small oscilla- 
tions in the density profile even in the presence of strong 
disorder, which is contrary to the repulsive case. In the 
inset, the amplitude of the ADW (Aadw), which is well 
defined close to the trap center, is shown for increasing 
disorder strength for two different values of \u\, u = —1 



and u = —2. The density distribution described here can 
be directly observed with the technique used by Foiling 
et ali^ They have observed the density distribution of a 
superfluid and a Mott insulator plateau for bosonic quan- 
tum gases in an optical lattice. We expect that the site 
occupation of fermionic systems will be directly detected 
in the near-future. 

The effects of disorder on the ADWs can be stud- 
ied through elastic light-scattering diffraction experi- 
ments?-^ by means of the Fraunhofer structure factor, 



S(k) 



1 

N] 



exp(— ikj)n 3 



The Fraunhofer structure factor has a maximum around 
the wave number k = kADW ~ ttL) with n being the 
average density in the bulk of the trap. Fig. [2] shows the 
structure factor for Nf = 30, N s = 100, u = — 2 and in 
the presence of a harmonic potential V^/i = 1.0 x 10~ 3 . 
Results for three cases of W are compared with the result 
for W = 0. In the system of W = 0, k ADW w 0.22, in 
which the Fraunhofer structure factor has a maximum. 
We can see how the ADW is suppressed by the disorder 
and thus S(k) is greatly suppressed accordingly. For W — 
0.5, S(k) becomes almost zero. 

We further analyze the effect of disorder on the tightly 
bound spin-singlet dimers at a strong coupling. In Fig.[3j 
we compare the density profile calculated from TEA with 
the one calculated by using DMRG. As it is clear, TFA 
gives us a very good overall shape for the strong attrac- 
tive interaction. The bound insulating region of tightly 
bound spin-singlet dimers at the center of the trap is 
quite unstable against disorder. This region shrinks 
quickly with an increasing disorder amplitude and van- 
ishes for a very small disorder value. The number of 
consecutive sites in the band-insulating range (iVdimer) 
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FIG. 3: (color online) DMRG results (open squares) for Nf — 
30, iV s = 200, Va/t = 1.0 x 10" 3 , and u = -20 in the clean 
limit are compared with TFA data (open circles). The thin 
solid lines aro just . a guide , for the ey e . 




FIG. 4: (color online) Site occupation Mi as a function of 
site position i for a system with Nf = 30 fermions in N s = 
200 lattice sites, and in the presence of a harmonic potential 
Va/t = 1.0 X 10 and strong attractive interaction of it = 
— 20. In the inset, the number of consecutive sites Ndimers 
at which |A/j — 2| < 10~ 5 is shown against the strength of 
the disorder W/t. The dashed line is the fitted exponential 
function: N dimera = -10.438 + 24.775e" 42 - 327(w ' /t)2 . 



is plotted as a function of the amplitude of disorder in 
the inset of Fig. QJ which shows that iVdimer = at 
W = 0.124. 

In Fig. O we illustrate the effect of disorder on the 
stiffness. The stiffness is defined as, S p = ((S^/SNf))^, 
which is the inverse of the global compressibility, and 
gives us precious information on the phases of the sys- 
tem. We would like to mention here, that the different 
quantum phases can be determined by time-of-flight ex- 
periments or by measuring noise correlations etc.£^ 

In the absence of disorder the ground state shows two 
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FIG. 5: (color online) Thermodynamic stiffness S p (in units of 
t) as a function of Nf for u = -20, V 2 /t = 10~ 3 and iV s = 200 
lattice sites. 



qualitatively different phases separated by one nonan- 
alyticity point: one phase characterized by the Luther- 
Emery liquid and the ADW, while the other with a band- 
insulator core of tightly bound spin-singlet dimers sur- 
rounded by Luther-Emery layers. The increase in stiff- 
ness is related to the incompressible nature of the insu- 
lating phase. 

The disorder has two main effects. Firstly it leads to an 
anomalous behavior in the stiffness at low density, which 
is similar to what has been found in the ID harmonically 
trapped Hubbard model of repulsive interactions,-- 5 - This 
phenomenon could be explained by the concept of den- 
sity percolation. The atoms tends to occupy the deepest 
valleys in the disorder landscape which makes the high- 
density regions disconnected in the low density limit. In 
other words, the system becomes more stiff. We con- 
clude that the low-density anomaly is robust against the 
nature of the interaction. Secondly, the disorder leads to 
the smearing of the nonanalytic point in the clean limit, 
signaling that the band-insulator core of tightly bound 
spin-singlet dimers is destroyed. 

In summary, on the basis of the Bethe-ansatz based lo- 
cal density-functional theory and the Thomas-Fermi ap- 
proximation, we have shown how disorder affects atomic 
density waves and spin-singlet dimers of attractively in- 
teracting Fermi gases in ID lattices under harmonic trap- 
ping potential. We have demonstrated that atomic den- 
sity waves are sensitive to the uniformly distributed dis- 
order. The tightly bound spin-singlet dimers in the cen- 
ter of the trap are quite unstable against the disorder. An 
anomalous increase of the stiffness persists at low density 
from quenching of percolation, which has also been found 
in repulsive systems. We have shown that the nonana- 
lyticity point in the stiffness signals the phase transition 
from the compressible liquid phase to the incompressible 
insulating phase with spin-singlet dimers in the center. It 
is found that very weak disorder leads to the smoothing of 
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the nonanalyticity point in the stiffness, thus suppresses 
the spin-singlet pairs. 
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